# encoding=utf-8




import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl


##############################################################################################

rf = 0.05                 # 无风险利率
mu = np.array([0.1,0.15,0.18]) # 预期收益率向量
sigma = np.array([0.1,0.12,0.15])
rou_ab, rou_ac, rou_bc = 0.2, 0.5, -0.2
# 协方差矩阵
C = np.array([[1,      rou_ab,  rou_ac],
              [rou_ab, 1,       rou_bc],
              [rou_ac, rou_bc,  1     ]])
V = np.outer(sigma,sigma)*C

# 前沿曲线绘图函数
def plot_efficient_frontier(rf=rf,mu=mu,V=V,is_subplot=False,show_E=False,strategy=None):
    f = mu-rf                 # 超额收益
    N = len(mu)
    e = np.ones(N)
    sigma = np.sqrt(np.diag(V))
    V_inverse = np.linalg.inv(V)
    # 等权投资组合
    w_E = e/N
    sigma2_E = w_E.dot(V).dot(w_E)
    sigma_E = np.sqrt(sigma2_E)
    f_E = w_E.dot(f)
    # 最小方差组合
    w_C = V_inverse.dot(e)/e.T.dot(V_inverse).dot(e)
    sigma2_C = w_C.dot(V).dot(w_C)
    sigma_C = np.sqrt(sigma2_C)
    f_C = w_C.dot(f)
    # 夏普组合
    w_Q = V_inverse.dot(f)/f.T.dot(V_inverse).dot(e)
    sigma2_Q = w_Q.dot(V).dot(w_Q)
    sigma_Q = np.sqrt(sigma2_Q)
    f_Q = w_Q.dot(f)
    # 有效前沿曲线
    f_P = np.linspace(f_C,1.5*np.max(mu)) # 只取上半部分抛物线
    get_w_P = lambda fp: (f_Q-fp)/(f_Q-f_C)*w_C + (fp-f_C)/(f_Q-f_C)*w_Q
    w_P = np.array([get_w_P(fp) for fp in f_P])
    k = (sigma2_Q-sigma2_C)/(f_Q-f_C)**2
    sigma2_P = sigma2_C + k*(f_P-f_C)**2
    sigma_P = np.sqrt(sigma2_P)
    r_P = rf + w_P.dot(f)
    # 策略组合
    if strategy is not None:
        w_S,f_S,sigma_S = strategy(pred_ret,V)
    # 资本市场曲线CML
    # s = np.linspace(0,0.2)
    # CML = rf + f_Q/sigma_Q * s
    # 有效前沿曲线绘制
    plt.plot(sigma_P,r_P)
#     plt.title('efficient frontier',fontsize=20)
    plt.title('efficient frontier',fontsize=10)
    plt.xlabel('$\sigma$',fontsize=20)
    plt.ylabel('return of Portfolio',fontsize=20)
    x_axis_low, x_axis_high = sigma_C-(sigma_P.max()-sigma_C)*0.1,sigma_P.max()
    plt.xlim(x_axis_low, x_axis_high)
    # 绘制风险资产
    plt.plot(sigma,mu,'bo')
    x_axis_length = x_axis_high-x_axis_low
    text_margin_space = x_axis_length*0.03
    for i in range(N):
        plt.text(sigma[i]+text_margin_space,mu[i]-text_margin_space,'asset '+str(i+1),fontsize=10)
    # 绘制等权投资组合
    if show_E:
        plt.plot(sigma_E, rf+f_E,'ko')
        plt.text(sigma_E+text_margin_space,rf+f_E-text_margin_space,'E',fontsize=15)
    # 绘制策略组合
    if strategy:
        plt.plot(sigma_S, rf+f_S,'ko')
        plt.text(sigma_S+text_margin_space,rf+f_S-text_margin_space,'S',fontsize=15)
    # 绘制最小方差组合与夏普组合（市场组合），以及无风险资产
    plt.plot(sigma_C, rf+f_C,'go')
    plt.text(sigma_C+text_margin_space,rf+f_C-text_margin_space,'C',fontsize=15)
    plt.plot(sigma_Q, rf+f_Q,'ro')
    plt.text(sigma_Q-text_margin_space,rf+f_Q+text_margin_space,'Q',fontsize=15)
#    plt.plot(0, rf,'yo')
#    plt.text(0-0.005,rf+0.008,'$r_{f}$',fontsize=15)
    # 绘制CML
    # plt.plot(s, CML, 'r-')
    # plt.text(0.15,0.28,'CML',fontsize=15)
    if not is_subplot:
        plt.show()


plt.figure(figsize=(15,9))
plot_efficient_frontier(show_E=True)


##############################################################################################
seed = 555
np.random.seed(seed)
# 真实收益率均值，用以生成数据，是不可观测的
true_mu = np.array([0.1,0.15,0.18])
true_sigma = np.array([0.1,0.12,0.15])
rou_ab, rou_ac, rou_bc = 0.25, 0.4, -0.3
true_C = np.array([ [1,      rou_ab,    rou_ac ],
                    [rou_ab,    1,      rou_bc],
                    [rou_ac,    rou_bc,   1   ]])
N = len(true_mu)
# 真实协方差矩阵，用以生成数据，是不可观测的
true_cov = np.outer(true_sigma,true_sigma)*true_C
def asset_generator(mean, cov, size=1):
    return np.random.multivariate_normal(mean, cov, size)
# 历史资产收益率历史
size = 20
asset_ret_obs = asset_generator(true_mu, true_cov, size=size)
N_obs = len(asset_ret_obs)
# 起初价格，可视化用
p0 = np.array([1,2,3])
# 还原历史价格序列
asset_price_obs = p0*np.exp(np.log(asset_ret_obs+1).cumsum(axis=0))
# 我们的收益预测模型
def mysterious_predict_model(asset_ret_obs,window=0):
    pred = asset_ret_obs[-window:,:].mean(axis=0)
    return pred
pred_ret = mysterious_predict_model(asset_ret_obs)
pred_price = asset_price_obs[-1,:]*(1+pred_ret)
estimate_cov = np.cov(asset_ret_obs.T) # 我们的协方差估计，不止这一种，这里只是示例
estimate_sigma = np.sqrt(np.diag(estimate_cov))
# 价格（－95%，95%）置信区间
pred_L, pred_H = asset_price_obs[-1,:]*(1+pred_ret-2*estimate_sigma), asset_price_obs[-1,:]*(1+pred_ret+2*estimate_sigma)

color_cycle = ['r', 'g', 'b']
# mpl.rcParams['axes.color_cycle'] = color_cycle


plt.figure(figsize=(15,5))
plt.subplot(2,2,1)
#for i in range(N):
plot1 = plt.plot(asset_ret_obs)
pred_pos=len(asset_ret_obs)
# 更新预测收益
for i in range(N):
    plt.plot([pred_pos-1,pred_pos],[asset_ret_obs[-1,i],pred_ret[i]],color_cycle[i]+'--')
plt.title('asset return series',fontsize=20)
plt.xlabel('time',fontsize=10)
plt.ylabel('return',fontsize=20)
plt.xlim(0,N_obs+1)
plt.legend(plot1,['asset1', 'asset2', 'asset3'],loc='upper left')
plt.tight_layout()

plt.subplot(2,2,2)
plot2 = plt.plot(asset_price_obs)
# 更新预测价格
for i in range(N):
    plt.plot([pred_pos-1,pred_pos],[asset_price_obs[-1,i],pred_price[i]],color_cycle[i]+'--')
    plt.plot([pred_pos,pred_pos],[pred_L[i],pred_H[i]],color_cycle[i]+'-.') # 绘制置信区间
plt.title('asset price series',fontsize=20)
plt.xlabel('time',fontsize=10)
plt.ylabel('price',fontsize=20)
plt.xlim(0,N_obs+1)
plt.legend(plot2,['asset1', 'asset2', 'asset3'],loc='upper left')
plt.tight_layout()

def our_strategy(pred_ret,V,rf=rf):     
    # 全额等额买入排名前50%的资产
    threshold = int(N/2)
    buy_list = (pred_ret.argsort()>=threshold)*1
    w_S = buy_list/buy_list.sum()
    sigma2_S = w_S.dot(V).dot(w_S)
    sigma_S = np.sqrt(sigma2_S)
    f_S = w_S.dot(pred_ret-rf)
    return w_S,f_S,sigma_S

plt.subplot(2,2,3)
# 后面会用到
mu_estimator = pred_ret
cov_estimator = estimate_cov
plot_efficient_frontier(rf=rf,mu=pred_ret,V=estimate_cov, is_subplot=True,strategy=our_strategy)

initial_steps = 20
simulated_steps = 100
size = initial_steps+simulated_steps
loops = 200

SR_S_list = []
SR_Q_list = []
for _ in range(loops):
    simulated_return = asset_generator(true_mu, true_cov, size=size)

    actial_rp_S_history = []
    actial_rp_Q_history = []
    for current_step in range(initial_steps,simulated_steps):
        curr_ret_obs = simulated_return[:current_step,:] # 当前的可观测数据集
        # 利用固定窗口数据，不然我们的模型会越来越准确，这里希望我们的预测模型的预测能力是一致的
        pred_ret = mysterious_predict_model(curr_ret_obs,window=20)
        estimate_cov = np.cov(curr_ret_obs.T) 
        estimate_sigma = np.sqrt(np.diag(estimate_cov))

        f = pred_ret-rf 
        N = len(pred_ret)
        e = np.ones(N)
        V = estimate_cov
        V_inverse = np.linalg.inv(V)
        # 策略组合 , 全额等额买入排名前50%的资产
        w_S,f_S,sigma_S = our_strategy(pred_ret,V)
        # 最小方差组合
        w_C = V_inverse.dot(e)/e.T.dot(V_inverse).dot(e)
        sigma2_C = w_C.dot(V).dot(w_C)
        sigma_C = np.sqrt(sigma2_C)
        f_C = w_C.dot(f)
        # 夏普组合
        w_Q = V_inverse.dot(f)/f.T.dot(V_inverse).dot(f) * sigma2_C * f.T.dot(V_inverse).dot(f)/f_C
        sigma2_Q = w_Q.dot(V).dot(w_Q)
        sigma_Q = np.sqrt(sigma2_Q)
        f_Q = w_Q.dot(f)
        #计算实际损益
        actial_return = simulated_return[current_step+1,:]
        actial_rp_S = w_S.dot(actial_return)
        actial_rp_Q = w_Q.dot(actial_return)
        actial_rp_S_history.append(actial_rp_S)
        actial_rp_Q_history.append(actial_rp_Q)
    SR_S = (np.mean(actial_rp_S_history)-rf)/np.std(actial_rp_S_history)
    SR_Q = (np.mean(actial_rp_Q_history)-rf)/np.std(actial_rp_Q_history)
    SR_S_list.append(SR_S)
    SR_Q_list.append(SR_Q)
plt.subplot(2,2,4)
plt.boxplot(x=[SR_S_list,SR_Q_list],labels=['S','Q'])
plt.title('Sharp ratio between S and Q',fontsize=10)
plt.xlabel('portfolio',fontsize=20)
plt.ylabel('sharp ratio',fontsize=20)

plt.show()







##############################################################################################




seed = 999
np.random.seed(seed)
# 有效前沿后验分布
Q_stats_list = []
C_stats_list = []

plt.figure(figsize=(15,3))
for _ in range(40):
    # 利用我们之前的参数估计模拟资产收益和协方差矩阵
    asset_ret_obs = asset_generator(mu_estimator, cov_estimator, size=size)
    pred_ret = mysterious_predict_model(asset_ret_obs)
    estimate_cov = np.cov(asset_ret_obs.T) 
    estimate_sigma = np.sqrt(np.diag(estimate_cov))
    
    f = pred_ret-rf 
    N = len(pred_ret)
    e = np.ones(N)
    V = estimate_cov
    V_inverse = np.linalg.inv(V)
    
    # 最小方差组合
    w_C = V_inverse.dot(e)/e.T.dot(V_inverse).dot(e)
    sigma2_C = w_C.dot(V).dot(w_C)
    sigma_C = np.sqrt(sigma2_C)
    f_C = w_C.dot(f)
    # 夏普组合
    w_Q = V_inverse.dot(f)/f.T.dot(V_inverse).dot(f) * sigma2_C * f.T.dot(V_inverse).dot(f)/f_C
    sigma2_Q = w_Q.dot(V).dot(w_Q)
    sigma_Q = np.sqrt(sigma2_Q)
    f_Q = w_Q.dot(f)
    # 有效前沿曲线
    f_P = np.linspace(f_C,1.5*np.max(pred_ret)) # 只取上半部分抛物线
    get_w_P = lambda fp: (f_Q-fp)/(f_Q-f_C)*w_C + (fp-f_C)/(f_Q-f_C)*w_Q
    w_P = np.array([get_w_P(fp) for fp in f_P])
    k = (sigma2_Q-sigma2_C)/(f_Q-f_C)**2
    sigma2_P = sigma2_C + k*(f_P-f_C)**2
    sigma_P = np.sqrt(sigma2_P)
    r_P = rf + w_P.dot(f)
    # 缓存Q、C组合统计信息
    Q_stats_list.append((rf+f_Q,sigma_Q))
    C_stats_list.append((rf+f_C,sigma_C))
    # 有效前沿曲线绘制
    plt.plot(sigma_P,r_P,'b')
# 绘制Q、C样本
plt.plot(np.array(Q_stats_list)[:,1],np.array(Q_stats_list)[:,0],'ro')
plt.plot(np.array(C_stats_list)[:,1],np.array(C_stats_list)[:,0],'go')
plt.title('efficient frontier posterior distribution',fontsize=20)
plt.xlabel('$\sigma$',fontsize=20)
plt.ylabel('return of Portfolio',fontsize=15)
plt.xlim(0.05,0.25)
plt.figure(figsize=(15,9))
plt.show()




##############################################################################################


# # 随机选取样本
# start_date = '2014-01-01'
# end_date = '2017-12-30'
# instruments = ['601601.SHA', '600030.SHA', '601668.SHA', '600036.SHA', '601098.SHA']
# df = D.history_data(instruments,start_date,end_date, fields=['open','close'])

# def get_logret(vec):
#     return np.log(vec).diff()[1:]

# # 重抽样还原股票相关数据
# resample_freq = 'M'
# df2 = df.set_index('date')[['instrument','close']]\
#                 .groupby('instrument')\
#                 .apply(lambda df: df['close'].resample(resample_freq,how='last'))\
#                 .transpose()
# ret_df = df2.apply(lambda x: get_logret(x))


# stock_mu = ret_df[:-5].mean()*12
# stock_V = ret_df.cov()*12

# rf=0
# mu=stock_mu
# V=stock_V
# f = mu-rf  
# N = len(mu)
# e = np.ones(N)
# sigma = np.sqrt(np.diag(V))
# V_inverse = np.linalg.inv(V)
# # 最小方差组合
# w_C = V_inverse.dot(e)/e.T.dot(V_inverse).dot(e)
# sigma2_C = w_C.dot(V).dot(w_C)
# sigma_C = np.sqrt(sigma2_C)
# f_C = w_C.dot(f)
# # 夏普组合
# w_Q = V_inverse.dot(f)/(f.T.dot(V_inverse).dot(e))
# sigma2_Q = w_Q.dot(V).dot(w_Q)
# sigma_Q = np.sqrt(sigma2_Q)
# f_Q = w_Q.dot(f)
# # 有效前沿曲线
# f_P = np.linspace(f_C,1) # 只取上半部分抛物线

# get_w_P = lambda fp: (f_Q-fp)/(f_Q-f_C)*w_C + (fp-f_C)/(f_Q-f_C)*w_Q

# w_P = np.array([get_w_P(fp) for fp in f_P])
# k = (sigma2_Q-sigma2_C)/(f_Q-f_C)**2
# sigma2_P = sigma2_C + k*(f_P-f_C)**2
# sigma_P = np.sqrt(sigma2_P)
# r_P = rf + w_P.dot(f)

# plt.figure(figsize=(15,5))
# plt.subplot(1,2,1)
# # 有效前沿曲线绘制
# plt.plot(sigma_P,r_P)
# plt.title('efficient frontier',fontsize=20)
# plt.xlabel('$\sigma$',fontsize=20)
# plt.ylabel('return of Portfolio',fontsize=20)

# plt.plot(sigma,mu,'bo')
# for i in range(len(instruments)):
#     plt.text(sigma[i],mu[i],str(instruments[i][:-4]))
# plt.plot(sigma_Q, rf+f_Q,'ro')
# plt.plot(sigma_C, rf+f_C,'go')

# plt.subplot(1,2,2)
# plt.bar(range(len(w_Q)),w_Q)
# plt.title("weight vector",fontsize=20)
# plt.ylabel("weight",fontsize=20)
# plt.xlabel("stock",fontsize=20)
# plt.xticks(range(len(w_Q)),[x[:-4] for x in instruments])

# plt.show()




##############################################################################################


# Q_portfolio_wealth = 1
# E_portfolio_wealth = 1
# Q_portfolio_wealth_path = []
# E_portfolio_wealth_path = []

# L = len(ret_df)
# start_index=24 # 预留24个月数据
# rf=0

# w_Q = None
# prev_w_Q = None
# trans_cost_ratio = 0.01
# N = ret_df.shape[1]
# e = np.ones(N)
# w_E = np.ones(N)/N
# for i in range(start_index,L):  
#     ## 每次循环的逻辑相当于期初开仓或调仓、期末计算收益，其中会考虑组合偏移因素
#     ## 即10%的资金开仓，期末由于股票上涨50%，导致组合权重偏移5%（假设其他股票没变化）
#     curr_ret_df = ret_df.iloc[i,:]
#     hist_2_prev_ret_df = ret_df.iloc[:i,:]
#     mu = hist_2_prev_ret_df[:-5].mean()*12 ## 用MA(5)作为下个时间段的预测
#     V = hist_2_prev_ret_df.cov()*12
#     f = mu-rf  
#     sigma = np.sqrt(np.diag(V))
#     V_inverse = np.linalg.inv(V)
#     # 夏普组合
#     adjusted_w_Q = V_inverse.dot(f)/(f.T.dot(V_inverse).dot(e))
#     r_Q = adjusted_w_Q.dot(curr_ret_df)
#     r_E = w_E.dot(curr_ret_df)
    
#     # 假设交易费用为净仓位的1%，包含交易佣金、手续费和滑点
#     if i == start_index:
#         c_Q = trans_cost_ratio*np.abs(adjusted_w_Q).sum()
#         c_E = trans_cost_ratio*np.abs(w_E).sum()
#     else:
#         # 计算组合偏移权重
#         shifted_w_Q = prev_w_Q * (1+curr_ret_df)
#         shifted_w_Q /= shifted_w_Q.sum()
#         shifted_w_E = w_E * (1+curr_ret_df)
#         shifted_w_E /= shifted_w_E.sum()
#         c_Q = trans_cost_ratio*np.abs(adjusted_w_Q-shifted_w_Q).sum()
#         c_E = trans_cost_ratio*np.abs(w_E-shifted_w_E).sum()
#     prev_w_Q = adjusted_w_Q
#     Q_portfolio_wealth *= 1+r_Q-c_Q
#     E_portfolio_wealth *= 1+r_E-c_E
#     Q_portfolio_wealth_path.append(Q_portfolio_wealth)
#     E_portfolio_wealth_path.append(E_portfolio_wealth)
#     #print(Q_portfolio_wealth, r_Q,c_Q, adjusted_w_Q.sum(),E_portfolio_wealth, r_E,c_E)

# plt.figure(figsize=(15,5))
# p1,=plt.plot(Q_portfolio_wealth_path)
# p2,=plt.plot(E_portfolio_wealth_path)
# plt.title('wealth curve',fontsize=20)
# plt.xlabel('time step',fontsize=20)
# plt.ylabel('account wealth',fontsize=20)
# plt.legend([p1,p2],['sharp portfolio', 'equal weight portfolio'],loc='upper left')

# plt.show()



##############################################################################################

